// STL
#include <fstream>                            // std::ifstream, std::ofstream
#include <memory>                             // std::unique_ptr<>
#include <string>                             // std::string
#include <utility>                            // std::move()
#include <vector>                             // std::vector<>

// Boost
#include <boost/archive/text_iarchive.hpp>    // boost::archive::text_iarchive
#include <boost/archive/text_oarchive.hpp>    // boost::archive::text_oarchive
#include <boost/serialization/unique_ptr.hpp> // serialize std::unique_ptr<>

// stats++
#include "statsxx/machine_learning/Gaussian_process.hpp" // Gaussian process stuff
#include "statsxx/statistics/CovarianceFunction.hpp"     // Covariance functions
#include "statsxx/statistics/covariance_functions.hpp"   // squared exponential


//
// DESC: Creates a GP (or reads one in).
//
std::unique_ptr<Gaussian_process::GaussianProcess> create_GP(
                                                             const bool        create,
                                                             // -----
                                                             const int         ni,
                                                             const int         no,
                                                             // -----
                                                             const std::string covariance_function,
                                                             const double      _theta,
                                                             // -----
                                                             const std::string gp_file
                                                             )
{
    std::unique_ptr<Gaussian_process::GaussianProcess> gp;

    std::cout << "here 1" << std::endl;

    if( create )
    {
        std::unique_ptr<CovarianceFunction> cf;

        // TODO: NOTE: Refresh memory of theta initialization (is it for all input and +1 for output?).

        std::vector<double> theta(ni, _theta);
        theta.push_back(1.);

        if( covariance_function == "squared_exponential" )
        {
            cf = std::unique_ptr<SquaredExponential>(new SquaredExponential(theta));
        }
        // NOTE: Correct specification of covariance_function should have been checked in the input.

        gp = std::unique_ptr<Gaussian_process::GaussianProcess>(new Gaussian_process::GaussianProcess(cf));
/*
        // SAVE GP TO ARCHIVE
        {
            std::ofstream ofs(gp_file);

            boost::archive::text_oarchive oa(ofs);

            oa << gp;
        }
*/
    }
    else
    {
        std::cout << "here 2" << std::endl;

//        std::unique_ptr<CovarianceFunction> cf;
//        Gaussian_process::GaussianProcess gp2(cf);

/*
        // READ GP FROM ARCHIVE
        {
            std::ifstream ifs(gp_file);

            boost::archive::text_iarchive ia(ifs);

            ia >> gp;
        }
        */
    }

    std::cout << "here 3" << std::endl;

    return std::move(gp);
}
